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Molecular dynamics simulations are performed to study spatially heterogeneous dynamics in a 
model of viscous silica above and below the critical temperature of the mode coupling theory, Tmct- 
Specifically, we follow the evolution of the dynamic heterogeneity as the temperature dependence 
of the transport coefficients shows a crossover from non-Arrhenius to Arrhenius behavior when the 
melt is cooled. It is demonstrated that, on intermediate time scales, a small fraction of oxygen and 
silicon atoms are more mobile than expected from a Gaussian approximation. These highly mobile 
particles form transient clusters larger than that resulting from random statistics, indicating that 
dynamics are spatially heterogeneous. An analysis of the clusters reveals that the mean cluster size 
is maximum at times intermediate between ballistic and diffusive motion, and the maximum size 
increases with decreasing temperature. In particular, the growth of the clusters continues when the 
transport coefficients follow an Arrhenius law. These findings imply that the structural relaxation 
in silica cannot be understood as a statistical bond breaking process. Though the mean cluster 
sizes for silica are at the lower end of the spectrum of values reported in the literature, we find that 
spatially heterogeneous dynamics in strong and fragile glass formers are similar on a qualitative level. 
However, different from results for fragile liquids, we show that correlated particle motion along 
quasi one-dimensional, string-like paths is of little importance for the structural relaxation in this 
model of silica, suggesting that string-like motion is suppressed by the presence of covalent bonds. To 
study spatial correlations between highly immobile particles, we calculate a generalized susceptibility 
corresponding to the self part of a four-point time dependent density correlation function. We find 
that this generalized susceptibility is maximum on the time scale of the structural relaxation, where 
a strong increase of the peak height indicates a growing length of spatial correlations between highly 
immobile particles upon cooling. Characterizing the local structures of the most mobile and the 
most immobile particles, respectively, we show that high particle mobility is facilitated by, but not 
limited to, the vicinity of defects of the network structure. 

PACS numbers: 66.30.Dn 
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I. INTRODUCTION 



Applications of various experimental methods demon- 
strated that the structural relaxation of many super- 
cooled liquids exhibits dynamic heterogeneities, i.e., it is 
possible to select a subset of particles that rotate or trans- 
late much farther or shorter distances than an average 
particlei£i2i2i£i& However, the majority of experimental 
techniques do not provide information about the spa- 
tial arrangement of mobile and immobile particles, which 
plays a central role in several theories of the glass tran- 
sitioni^iS Recently, results of multidimensional nuclear 
magnetic resonance experiments showed that dynamics 
in polymersifi and in organic low-molecular weight com- 
poundsAi*i£ are spatially heterogeneous. Particles within 
a physical region of these viscous liquids show an en- 
hanced or diminished mobility compared to particles in 
a region a few nanometers away. Despite this progress, 
an experimental characterization of the time evolution of 
spatially heterogeneous dynamics (SHD) is still lacking. 
Further, detailed experimental studies of the dynamical 
behavior of strong liquids, like silica, are hampered by 
the high glass transition temperatures, T g , of these ma- 
terials. Nevertheless, based on the Arrhenius behavior 
observed at sufficiently low temperatures, it is often ar- 



gued that the dynamics in silica melts can be understood 
as a statistical bond breaking process. 

Molecular dynamics (MD) simulations of model glass- 
forming liquids provide direct access to spatial correla- 
tions of particle mobility. Due to the available com- 
puter power, present simulations are performed at tem- 
peratures near Tmct, the critical temperature of the 
mode-coupling theory (MCT) ji 3 . where the structural re- 
laxation typically occurs on the order of nanoseconds. 
In previous work, MD simulations were used to inves- 
tigate dynamic heterogeneity in soft-spherai^iiSii& and 
hard-sphere systems^ binary Lennard- Jones (BLJ) mix- 
tures, 18 . 19 ! 20 : 21 ! 22 . 23 ! 24 ! 25 ! 26 ! 27 the Dzugutov liquid,^?- 3 ^ 
polymer melto 32 ' 33 : 34 : 35 and molecular liquidsi 36 ' 37 These 
studies showed that highly mobile and highly immobile 
particles aggregate into clusters that are transient in na- 
ture. 

Specifically, highly mobile particles both in simple liq- 
uids, such as BLJ mixtures 2 ^ and the Dzugutov liquid, 30 
and in more complex systems, like polymer melts 3 ^ and 
water, 3 ' form clusters that are largest at times inter- 
mediate between ballistic and diffusive motion. The 
size of these clusters grows strongly when the liquid is 
cooled and a divergence at T^Tmct was proposed. 22 : 38 
Within the clusters smaller groups of particles move in 
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strings, i.e., particles follow each other along quasi one- 
dimensional paths. This correlated motion on interme- 
diate time scales, which is believed to facilitate struc- 
tural relaxation,- 8 appears to be universal as well, as it 
is found for a range of different model liquids iSLSS^SiS 
In addition, spatial correlations of highly immobile par- 
ticles were studied based on a four-point time dependent 
density correlation functiom 24 i 25 i 26 i 27 i 28 i 39 For BLJ liq- 
uids, it was observed that the corresponding generalized 
susceptibility X4(Ai) exhibits a maximum that increases 
with decreasing temperature, indicating that there is a 
growing length of spatial correlations between particles 
localized in a time interval At. 

All these studies were performed in temperature ranges 
where the evolution of the relaxation times shows pro- 
nounced deviations from Arrhcnius behavior. In a recent 
MD simulation study^i we demonstrated that SHD ex- 
ists not only in fragile liquids, but also in a model of 
silica, the paradigm of a strong liquid. 41 Specifically, the 
most mobile particles were found to form clusters larger 
than predicted from random statistics, where the mean 
cluster size is maximum at times intermediate between 
ballistic and diffusive motion, as was observed for other 
model liquids. Further, we showed that the dynamic het- 
erogeneities are short lived so that high particle mobil- 
ity spreads throughout the sample on the time scale of 
the structural relaxation. In doing so, the probability 
for a previously immobile particle to become mobile is 
enhanced in the vicinity of another mobile particle, sup- 
porting the concept of dynamic facilitation, which is the 
cornerstone of a microscopic model of viscous liquids pro- 
posed recently by Garrahan and Chandler. 8,9 Finally, our 
previous results suggest that, compared to other model 
liquids, string-like motion is less relevant for the struc- 
tural relaxation of silica. 

In view of these findings for fragile and strong liquids, 
one may ask how the transport coefficients of viscous liq- 
uids - in particular their temperature dependence - are 
related to the properties of SHD. Indeed, the decoupling 
of diffusivity and viscosity or structural relaxation time, 
indicating a breakdown of the classic Stokes- Einstein re- 
lation, has been rationalized in terms of SHDi£*2i Fur- 
thermore, a link between the diffusion coefficient D and 
the clusters of mobile particles was observed in simula- 
tion studies of water £L Specifically, it was reported that 
the mean cluster size is a measure of the mass of the coop- 
eratively rearranging regions central to the Adam-Gibbs 
theory 7 which in turn allows one to describe the temper- 
ature dependence of D. This suggests that the behavior 
of the transport coefficients, which describe the long-time 
liquid dynamics, is determined by SHD at intermediate 
times. On the other hand, Garrahan and Chandler^ 
emphasize that spatiotemporal characteristics of mobil- 
ity propagation are important to rationalize differences 
between various glass formers. In particular, they argue 
that mobility propagation carries a direction, the persis- 
tence length of which is larger for fragile than for strong 
glass formers. 



Here, we continue our MD simulation study of SHD 
in the BKS model of viscous silica^ which is commonly 
used to reproduce structural and dynamical properties 
of this strong liquidi 43 i 44 i 45 i 46 i 47 i 48 i 49 i 50 i 51 In this way, we 
wish to obtain further insights into the relaxation dy- 
namics of one of the most important glass formers, which 
are difficult to extract from experimental studies. It is 
well established that the transport coefficients of BKS sil- 
ica show a crossover from a non-Arrhenius to an Arrhe- 
nius temperature dependence upon cooling 44 i 45 i 46 i 47 i 48 
in agreement with experimental data. 52 Saika-Voivod et 
al4& related such behavior to a "fragile-to-strong transi- 
tion" . A main goal of the present work is to quantify the 
evolution of SHD during this crossover. Moreover, by 
means of quantitative comparison with literature data, 
we seek to ascertain differences and similarities of SHD 
among fragile and strong liquids, respectively. We study 
spatial correlations both between highly mobile particles 
and between highly immobile particles. Finally, we ana- 
lyze the relation between properties of the local instan- 
taneous liquid structure and the particle mobility, to in- 
vestigate whether there is a structural origin of SHD in 
silica. Indications of a relationship between local struc- 
ture and local dynamics have been demonstrated in the 
DzugutoMS 9 ^ 3 . and BLJ liquids^ 

II. MODEL AND SIMULATION 

Previous MD simulation studieo 43 i 44 i 45 i 46 i 47 i 48 i 49 i 50 i 51 
have shown that the BKS potential 4 ^ is well suited to 
reproduce many structural and dynamical properties of 
amorphous silica. The BKS potential energy is given by 

4><xp{r) = — h A a[} exp(-B al3 r) — (1) 

where r is the distance between two atoms of types a and 
(3 (a, j3 € {Si, O}) and e denotes the elementary charge. 
The partial charges q a together with the potential pa- 
rameters Aa/3, B a j3 and C a p can be found elsewhere 

Here, we follow simulations of Horbach and Kob 44,45 
and perform computations in the NVE ensemble where 
N = 8016 and p = 2.37 g/cm 3 corresponding to a box 
length L = 48.4A. This system size is sufficiently large 
so as to avoid significant finite size effects at the studied 
rp 50,51 rp-^g e q ua tions of motion are integrated using the 
velocity Verlet algorithm with a time step of 1.6 fs. In 
doing so, we truncate the non-Coulombic part of the po- 
tential at a cutoff radius of 5.5 A. The Coulombic part is 
calculated via Ewald summation 5 * where we apply a pa- 
rameter a = 0.265 together with a cutoff in Fourier space 
k c = 27r/Lv / 5T. As was discussed in previous workj 4 ^ 5 - 
the use of this value of k c leads to a constant shift of 
the total potential energy and, hence, does not affect 
the forces. We consider values of T in a range between 
3030 K and 5250 K, where the system is equilibrated for 
times longer than the structural relaxation time prior to 
data accumulation. 
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Horbach and Kob 44,45 showed that BKS silica consists 
of a network of well defined silicate tetrahedra. Further, 
they analyzed the temperature dependences of the diffu- 
sion coefficient D and the time constant t of the struc- 
tural relaxation. At relatively high T, these quantities 
were found to vary according to a power-law 

D, t- 1 ex {T-T MCT )\ (2) 

as predicted by MCTpi^ For the critical temperature, 
Horbach and Kob^ii 5 . reported Tmct — 3330 K, in good 
agreement with Tmct — 3221 K determined from exper- 
imental viscosity data of silica^ The critical exponents 
extracted from D{T) and r _1 (T) were found to be7«2.1 
and 7 w 2.4, respectively— 4 s . At T below a dynamical 
crossover in the vicinity of Tmct , Arrhcnius laws 

D, t- 1 cx exp [-E a /{k B T)\ (3) 

were observed Activation energies E a =4.7eV and 
E a = 5.2 eV for the oxygen and silicon atoms, respec- 
tively, were obtained from D(T). The temperature de- 
pendence r _1 (T), extracted from the incoherent interme- 
diate scattering functions for the oxygen atoms, yielded 
activation energies between 5.0 eV and 5.5 eV, depend- 
ing on the wave vector q. All these values are similar 
to E a — 4.7 eV for oxygen and E a = 6.0 eV for silicon 
observed in experimental studies near T g — 

Our analysis confirms these findings of Horbach and 
Kob— ^ Specifically, our data are consistent with 
Tmct = 3330 K and we obtain critical exponents 7«2.0 
and 7 sa 2.2 from D(T) and r~ 1 (T), respectively. Fur- 
thermore, we observe comparable activation energies, see 
section 1111 Al 

III. RESULTS 

A. Properties of the bulk 

To mark the various time regimes of the structural re- 
laxation in BKS silica, we first discuss the incoherent 
intermediate scattering function 

F s (q,At) = (cos{q[n{t + At) - r t {t )}}). (4) 

Here, (to) is the position of particle i at time to and 
the brackets (. . .) denote the ensemble average. Figure^] 
shows F s (q, At) for the oxygen and silicon atoms, respec- 
tively, where we use an absolute value of the wave vec- 
tor, q — 1.7A . This value of q corresponds to the first 
sharp diffraction peak of the static structure factor— and, 
hence, dynamics on the length scale of the distance of the 
silicate tetrahedra is probed. We see in Fig. that the 
scattering functions for both atomic species are compa- 
rable. Typical of viscous liquids, F s (q,At) shows a two- 
step decay. While the short-time decay can be ascribed 
to ballistic motion, the non-exponential long-time decay 
results from the structural relaxation. As a consequence 
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FIG. 1: Incoherent intermediate scattering functions, F s (q = 
1.7A ; At), for the oxygen and silicon atoms, respectively. 
The temperatures are, from left to right: 5250 K, 4730 K, 
4330 K, 4120 K, 3870 K, 3710 K, 3520 K, 3420 K, 3330 K, 
3230 K, 3140 K and 3030 K. 



of the cage effect^ which describes that the particles are 
temporarily trapped in a cage formed by their neighbors, 
a plateau regime develops between the ballistic and dif- 
fusive regimes upon cooling. At the crossover from the 
ballistic to the plateau regime, the curves show oscilla- 
tions that are damped out for longer times. Horbach et 
al. 50,51 found a similar behavior and ascribed it to the 
Boson peak. However, the origin of this feature is still a 
matter of debate— *5£ 

For an analysis of the structural relaxation, we fit the 
long-time decay of F s (q,At) to a Kohlrausch- Williams- 
Watts (KWW) function, A exp[-(t/r) /3 ]. In Fig. H it 
is evident that the temperature dependence of the mean 
time constant, given by (r) = (t / (3)T (1 / f3) , is similar 
for both atomic species. At T < 3300 K, the data are 
well described by an Arrhenius law where we find activa- 
tion energies E a — 5.0eV and E a = 5.1eV for the oxygen 
and silicon atoms, respectively. In contrast, there are 
deviations from an Arrhenius law at higher T, reflect- 
ing the crossover observed for the temperature depen- 
dence of the transport coefficientsi^^ 4 *^ In particular, 
(t (T)) can be fitted by Eq. H with Tmct = 3330 K at 
T > 3900 K. The non-exponentiality of the structural re- 
laxation can be quantified by the stretching parameter (3 
of the KWW function. We observe (3 > 0.8 for all studied 
T and for both atomic species, and, hence, the devia- 
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FIG. 3: Self part of the van Hove correlation function, 
G s (r,At), for the oxygen atoms in BKS silica at T = 3230K. 
The time At = 31.8 ps compares to the transition between 
the plateau and diffusive regimes. Dashed line: Gaussian 
approximation Go(r, At) for At = 31.8 ps, cf. Eq. |S| The 
arrows indicate the distances r c for which the integral (j> = 
f°° 47rr 2 G a (r,At = 31.8ps)dr equals 3%, 5% and 7%, respec- 
tively. The data to the right of the arrows result from highly 
mobile particles studied in section UlI Bl For comparison, the 
oxygen-oxygen interparticle distance roo equals 2.6 A, see 
Fig. [TT] 



FIG. 2: Various time constants characterizing the dynamics 
of the oxygen and silicon atoms in BKS silica, see text for 
details. The data (r (T)) are fit by a power-law at T > 3900 K 
(dashed lines, Eq. [5] with Tmct — 3330 K, 7 = 2.2 for oxygen 
and silicon) and by an Arrhenius law at T < 3300 K (solid 
lines, Eq. |3] -B a = 5.0eV for oxygen, _E a = 5.1eV for silicon). 

tions from an exponential behavior are small. A closer 
inspection reveals that starting from /jwl at high T the 
stretching parameter first decreases, and then becomes 
constant within statistical error at T < 3800 K (oxygen: 
/3 = 0.83±0.02, silicon: /3 = 0.85±0.02). 

Next, we study the self part of the van Hove correlation 
function 

G s (f,At) = (^[^(to + AtJ-^to) -»*]>. (5) 

The quantity 47rr 2 G s (r, At) measures the probability 
that a particle moves a distance r in a time interval At. 
Figure |21 depicts this probability for the oxygen atoms at 
T = 3230K. It is instructive to compare G s (r, At) with 
the Gaussian approximation^S* 2 ^^ 

Go(r ' At) = i^wm) 2 exp {-m*®) (6) 

where rf(At) = \n (t + At) - n(to)\ 2 . While G s (r, At) = 
Go (r, At) for short and long times in the structural relax- 
ation process, significant deviations are obvious at inter- 
mediate times. They are demonstrated for At = 31.8 ps 
in Fig. |3| For this time interval, Airr 2 G s (r, At) shows 
a pronounced tail, indicating that some particles move 



much farther than expected from a Gaussian approxima- 
tion. Such behavior is well known from previous studies 
of fragile viscous liquidsi 20 i 22 i 30 i 34 i 35 i 36 i 37 In the present 
system, we observe this behavior also at T below the 
fragile-to-strong crossover, consistent with previous find- 
ings^ This result gives a first indication of a similarity 
in SHD between fragile and strong liquids. 

The deviations from Gaussian behavior can be quan- 
tified by the non-Gaussian paramete r! 22 ' 35 

3 (rj(At)) 
2( ) = 5(r 2 (At)) 2_L (?) 

In Fig. 01 we see that <X2(At) exhibits a maximum at 
intermediate times. For both atomic species, the posi- 
tion of the maximum, t a i , shifts to longer times and the 
maximum value, o?2=a<i(t a <x), increases with decreasing 
T. The insets of Fig. 0] demonstrate that the increase 
of a™ can be described by an exponential growth with 
l/T. For the silicon atoms, the temperature dependence 
may be weaker at T <Tmct, but statistics are too poor 
to draw unambiguous conclusions. In any case, the data 
show that the crossover observed for the transport co- 
efficients is not of great importance for the temperature 
dependence of a™. At any given T, the non-Gaussian 
parameter is larger for oxygen than for silicon, implying 
that dynamic heterogeneity is more pronounced for oxy- 
gen. The temperature dependence of t a 2 is displayed in 
Fig. |21 Upon cooling, t a2 decouples from (r) so that 
the time constants differ by about one order of mag- 
nitude at T — 3030 K. Comparison with Fig. ^ shows 
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FIG. 4: Non-Gaussian parameter ct2(At) for the oxygen 
and silicon atoms at various temperatures (5250 K, 4330 K, 
3870 K, 3520 K, 3330 K, 3140 K, 3030 K). The insets show the 
temperature dependence of the maximum values, a™. The 
solid lines are interpolations of the data at T > 3600 K with 
an Arrhenius law. The limited fitting range was chosen to 
determine whether a™ shows a significantly different temper- 
ature dependence at low temperatures. 



that the deviations from Gaussian behavior are largest 
in the late-/?/ early-a relaxation regime, consistent with 
previous results for various models of supercooled liq- 
uid: . 20 - 22 - 30 - 34 - 35 - 36 - 37 



B. Properties of highly mobile particles 

In this section, we investigate the properties of highly 
mobile particles on intermediate time scales. Since we 
are interested in a detailed comparison with results for 
non-networked glass forming liquids, our analysis focuses 
on quantities studied in previous wor ] tI 22i22iMi2L 3 £ There, 
fractions of mobile particles <fr = 5 — 7% were considered. 
These fractions were obtained as the fractions of particles 
in the tail of G s (r,t a 2). Note that these highly mobile 
particles show displacements of at least one interparticle 
distance, while the rest of the particles are trapped in 
their local cages, see Fig. |3J Here, we compare results 
for <j> = 3%, (f> = 5% and <j> = 7%, where the most mobile 
particles in a time interval At are identified based on the 
scalar particle displacements within this time window. 
We analyze the behavior of the oxygen and silicon atoms 




1000 



1 10 
At [ps] 



1000 



FIG. 5: Normalized number-averaged mean cluster size, 
Sat (At), for the oxygen and silicon atoms at various tem- 
peratures (5250 K, 4330 K, 3870 K, 3520 K, 3330 K, 3140 K, 
3030 K) and a fraction = 5%. The insets show the probabil- 
ity distribution of the cluster size, Ps(n; At), for T = 3030K 
and T = 4330 K, and times when the mean cluster size is a 
maximum (At = ts)- 

separately, since the atomic species exhibit somewhat dif- 
ferent properties, as shown below. 



1. Clusters 

To ascertain the spatially heterogeneous nature of dy- 
namics in BKS silica we first show that highly mobile 
particles form clusters larger than expected from random 
statistics. Following previous studies; 22 i 30 i 34 i 37 i 38 i 40 we 
define a cluster as a group of highly mobile particles that 
reside in the first neighbor shells of each other. For both 
atomic species, we define the neighbor shell based on the 
first minimum of the respective pair correlation function. 
A statistical analysis of the clusters is possible when we 
determine the probability distribution Ps(n; At) of find- 
ing a cluster of size n for a time interval At. From this 
distribution, the number-averaged and weight-averaged 
mean cluster size can be calculated according to 



and 



S N {At) = nP s (n;At) 

n 

n 2 P s (n;Ai) 



Sw(At) 



J2 n nP s (n-At) ' 



(8) 



(9) 
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respectively, where J2n Ps( n ) = 1. 

In the following analysis of the clusters and, later, of 
the strings, we focus on the number-averaged data for 
the fraction <fi = 5%. Here, the number average is con- 
sidered due to the smaller statistical error of this quan- 
tity. However, we carefully checked that our conclusions 
depend neither on the average nor the fraction consid- 
ered. For several examples, this is demonstrated by in- 
cluding results for = 3% and = 7%. Some findings for 
the weight-averaged data were shown in previous work^S 
They are discussed in section lTVl when we compare SHD 
in various models of viscous liquids on a quantitative 
level. Finally, we determined that an analysis of clus- 
ters consisting of both mobile oxygen and silicon atoms 
does not yield new insights. 

In Fig. El we display the temperature dependence of 
the normalized mean cluster size S'jv(Ai) = Sjv(Ai)/S^-, 
where = 1.21 is the mean cluster size resulting when 
5% of the particles arc randomly chosen to construct clus- 
ters. Thus, S'Ar(At) measures exclusively effects due to 
SHD. Inspecting the oxygen data, we see that SV(Ai) 
exhibits a peak, indicating the existence of clusters that 
are transient in nature. This peak grows and shifts to 
longer times upon cooling. In addition, there is a shoul- 
der near Ai«0.3ps. These findings confirm our previ- 
ous results for the corresponding weight-averaged data, 
SV(A£)4°. For silicon, the shoulder becomes a separate 
maximum and Sjv(Ai) shows two peaks at low T. Since 
the position and height of the primary maximum de- 
pend more strongly on temperature than those of the 
secondary maximum, both peaks merge at high T. Com- 
paring the data for both atomic species, we observe that 
the clusters are larger for oxygen, consistent with a larger 
value of a™ for this species. 

A shoulder in SV(Ai) at the crossover from the bal- 
listic to the plateau regime was also observed for wa- 
ter^ whereas such an increase is absent in simple liq- 
u j c j 9 22 i 30 J 38 anc j a p ]y mer melti^i In the case of water j2I 
this phenomenon was ascribed to "strong correlations in 
the vibrational motion of first-neighbor molecules, owing 
to the presence of hydrogen bonds" . Similarly, we suggest 
that the network structure of silica enables a collective 
vibrational mode, which becomes more relevant when T 
is decreased. This conclusion is corroborated by the os- 
cillatory behavior of F s (q, At) in the corresponding time 
regime, see Fig. ^ Provided these findings for silica can 
be attributed to the Boson peak, as proposed by Hor- 
bach et aL^2*Si our findings imply that this phenomenon 
results from a local, collective vibration. 

SHD related to the structural relaxation manifests it- 
self in the primary maximum of S'jv(At). For a quantita- 
tive analysis, we extract the peak times and peak heights. 
In Fig. [21 we see that the mean cluster size is maximum 
at times ts -C (r), where the ratio (r)/ts is nearly in- 
dependent of T. For the silicon atoms, such analysis is 
not possible at sufficiently high T, because the primary 
maximum is submerged by the secondary maximum so 
that the temperature-independent position of the latter 
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FIG. 6: Maximum of the normalized mean cluster size, S N , 
together with the maximum of the mean string length, L™, for 
the oxygen and silicon atoms as a function of the reciprocal 
temperature. The various fractions of mobile particles, 0, 
used in the analysis are indicated. The solid lines are guides 
to the eye. 



results in a constant value for ts- Despite slight devi- 
ations, the temperature dependence of ts for different 
values of is comparable. Specifically, ts increases with 
increasing percentage, but the data for = 3% and = 7% 
differ by less than a factor of 2 at all studied T. In Fig. El 
we show the temperature dependence of the peak height 
S N = Siy(ts) for all considered fractions 0. Inspecting 
the data for the lower T, we observe that S 
nearly exponentially with increasing 1 /T 
a temperature range where the relaxation time (r) varies 
according to an Arrhenius law, the cluster size strongly 
increases. A more detailed analysis together with a dis- 
cussion of the results will be presented in section llH CI 

Finally, we investigate the probability distribution of 
the cluster size, Ps{n; At). For simple Iiquids224^2i2£ and 
a polymer melt^ this distribution exhibits a power-law 
behavior when At = ts and T « Tmgt- I n the insets 
of Fig. El we show Ps(n;ts) for the oxygen and silicon 
atoms at representative values of T. At T = 3030 K, 
clusters containing as many as 40 oxygen atoms and 20 
silicon atoms are observed. However, a power-law is not 
found in the studied temperature range, indicating that 
this feature of SHD in simple liquids cannot be general- 
ized to the case of the network-former BKS silica. The 



N grows 
Thus, even in 
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data can be satisfactorily fit by a power-law multiplied by 
an exponential cutoff. Such behavior is expected when a 
percolation transition is approached and, hence, we can- 
not exclude that such a transition occurs at much lower 
T. 



2. Strings 

It has been demonstrate d 21 - 22 -??'?^ that dynamics in 
several fragile liquids are facilitated by string-like mo- 
tion, i.e., groups of particles follow each other along one- 
dimensional paths. However, our previous work^ii sug- 
gests that this type of motion is less important for silica. 
To study the relevance of string-like motion in more de- 
tail, we follow Donati et al~ and construct strings by 
connecting any two particles i and j of the same atomic 
species if 

mm[\n{t a )- fj (t + At) | , | n (t + At) - P s (t ) | ] < 5. 

Similar to the values used in previous we 
set 5 to ~ 55% of the respective interatomic distance, 
resulting in 5 = 1.4 A and S = 1.7 A for oxygen and sili- 
con, respectively. Then, the above condition implies that 
one particle has moved and another particle has occu- 
pied its position. We checked that our conclusions are 
not affected when S is varied in a meaningful range. 

For an analysis of string- like motion, we determine the 
probability distribution Pl(1; At) of finding a string of 
length I in a time interval At and calculate the number- 
averaged mean string length L^(At) in analogy with 
Eq.HO Note that we include "strings" containing only one 
atom. In Fig. we display the temperature dependence 
of Lpf(At). Though the strings for both atomic species 
grow and shrink in time, string-like motion is very differ- 
ent for the oxygen and silicon atoms. We see that, at any 
given T, Ljv(At) is maximum at a much later time for 
silicon than for oxygen. Further, the mean string length 
is significantly smaller for the former than for the latter 
atomic species. In particular, for the silicon atoms, the 
small values L^(At) ~ 1 imply very limited string- like 
motion. 

The different behavior of the oxygen and silicon atoms 
can be quantified by the peak times t l and peak heights 
Ljy = Liv(tL)- The former are included in Fig. |21 For oxy- 
gen, we find that the mean string length and the mean 
cluster size are maximum at similar times t$~tL<^ (r), 
as was observed for simple liquidsS 2 ^ and a polymer 
melt^ In contrast, these mean sizes peak at very dif- 
ferent times for silicon where the few short strings are 
largest in the a-relaxation regime, i.e., ts <C tj, ~ (r). 
This finding clearly shows that the string-like motion 
known for fragile liquids is suppressed for the silicon 
atoms. In Fig. El we show the temperature dependence of 
the peak heights for all studied fractions 0. As was 
observed for the clusters, the mean string size for the 
oxygen atoms increases nearly exponentially with 1 /T at 
the lower T. For the silicon atoms, the strings grow more 
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FIG. 7: Number-averaged mean string length, Ljv(Ai), for 
the oxygen and silicon atoms at various temperatures (5250 K, 
4330 K, 3870 K, 3520 K, 3330 K, 3140 K, 3030 K), where = 
5%. The insets show the probability distribution of the string 
length, P(l;At), for T = 3030K and T = 4330K, and times 
when the mean string size is a maximum (At=ti). 

slowly, but large scattering of the data, in particular at 
low T, hampers a quantification of the temperature de- 
pendence. 

Examples of the probability distribution of the string 
length, Pl(1, At), are shown in the insets of Fig.[7| For all 
considered T and both atomic species, Pl(1, At) decays 
exponentially when At =tf, , as was observed in previous 
studies of string-like motioni^L 2 ^ ^ This finding sug- 
gests that an exponential distribution P^ti^l) is a fea- 
ture inherent to strings constructed in the way described 
above. In view of the exponential distribution, an anal- 
ogy to the equilibrium polymerization of linear polymers 
was proposed in previous work^i^ 

Finally, we analyze the relevance of string-like motion 
for the relaxation of the most mobile particles in silica. 
Following previous work^S we quantify the importance 
by the fraction of mobile particles, /(At), that are in- 
volved in non-trivial strings, i.e., strings with I > 3. In 
Fig. |S1 we see that string-like motion is of very limited 
relevance for silicon, whereas it becomes increasingly im- 
portant for oxygen upon cooling. Hence, this dynamical 
pattern may be an important channel of relaxation for 
oxygen near T g . A more detailed analysis reveals that 
/(At) is maximum at times At = tf « tj, for both the 
oxygen and the silicon atoms, see Fig. [21 Inspecting the 
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FIG. 8: Fraction f(At) of highly mobile oxygen and silicon 
atoms, respectively, that are involved in non-trivial string-like 
motion (Z>3). We display results for 5250 K, 4330 K, 3870 K, 
3520 K, 3330 K, 3140 K and 3030 K. The insets show the max- 
imum value, f m , as a function of the reciprocal temperature. 
The values of 4> used in the analysis are indicated. 



maximum amplitude f m in the insets of Fig. [8] it is evi- 
dent that our conclusions do not depend on the value of 



C. Relation to the Adam-Gibbs theory 

The Adam-Gibbs (AG) relation 7 . D oc exp[-A/(TS c )], 
which links the diffusion coefficient with the configura- 
tional entropy S c , was found to hold for a BLJ liquid, 58 
water— and the present case of BKS silica— The AG 
theory further proposes a relation between S c and the 
characteristic mass of cooperatively rearranging regions 
(CRR). However, the CRR are not precisely defined in 
the theory and their nature is still elusive. In recent work 
on water, Giovambattista et o/^i observed that 



(10) 



suggesting that, first, the (non-normalized) mean cluster 
size is a measure of the mass of the CRR and, second, a 
cluster of size one does not correspond to a CRR. Fur- 
thermore, they confirmed along the considered isochore 
the expectation 
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FIG. 9: Sff—1 and L™— 1 as a function of the logarithm of the 
diffusion coefficient D. The dashed lines are guides to the eye. 
The solid line is a linear interpolation, which will be expected 
if clusters or strings are related to the cooperatively rearrang- 
ing regions of the Adam-Gibbs theoryp— *— In the insets, we 
depict the temperature dependence of Xn — 1 = Sm— 1, L™— 1 
together with Arrhenius fits (solid lines). 



D<xexp[-A*(S% -l)/T] 



(11) 



resulting from the AG relation together with Eq. 

Here, we analyze whether Eq.^^is valid for BKS silica. 
Since the AG relation holds for this model in the stud- 
ied temperature ranged the validity (failure) of Eq. El 
implies the validity (failure) of Eq. ^| and, hence, we 
can determine whether a relation between the clusters 
of mobile particles and the CRR exists, as was reported 
for water— To this end, we extract the diffusion coef- 
ficient from the long-time behavior of the mean square 
displacement and plot (S^—\)/T as a function of log D 
in Fig. |5J Evidently, a linear relation is observed for nei- 
ther oxygen nor silicon, indicating a failure of Eq. 1111 
and, thus, also Eq. ^3 Hence, for silica, the mean clus- 
ter size is not a measure of the characteristic mass of 
the CRR. We also find that such relation is not valid for 
the weight-averaged data Sy^ or for clusters consisting 
of both oxygen and silicon atoms. On the other hand, 
it has been argued^Siffi that the strings and not the 
clusters are the elementary units of cooperative motion 
and, hence, related to the CRR. Therefore, we include 
{U$-l)/T as a function of logD in Fig.^J While no lin- 
ear relation is found for the oxygen atoms, such behavior 
cannot be ruled out for the silicon atoms. However, the 
variation of is too small to draw any conclusion in 
the latter case, where string-like motion is peculiar and 
largely insignificant anyway. 

In the insets of Fig. EI we show the temperature depen- 
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FIG. 10: Generalized susceptibility Xsa(At) for the oxygen 
and silicon atoms. The temperatures are 5250 K, 4330 K, 
3870 K, 3330 K, 3140 K and 3030 K. 



dence of S^ — l and L^ — 1. For both atomic species, we 
observe that an exponential growth with 1/T fits the data 
well over the entire temperature range. This behavior is 
confirmed by our analysis for the fractions <fi = 3% and 
= 7% (not shown), indicating that, in a strong glass for- 
mer, SHD becomes increasingly important with decreas- 
ing T. In particular, we find no evidence for substantial 
effects due to the "fragile-to-strong crossover" observed 
for the transport coefficients. Provided the exponential 
growth continues upon further cooling, the mean cluster 
and mean string size are finite at all finite T . In con- 
trast, a power-law temperature dependence of the mean 
cluster size was found for a BLJ liquid ) 22 i 38 suggesting 
a percolation transition of clusters of mobile particles in 
the vicinity of Tmct- 



D. Behavior of four-point spatiotemporal density 
fluctuations 

Previous work^iiS^iSSiSLiS showed that two-point, two- 
time fourth order density correlation functions are well 
suited to study SHD. The general theoretical framework 
is described in detail in the literaturei^S^ In brief, it has 
been proven useful to define an "order parameter" Q(At) 
that compares the configurations of a liquid with density 



p(f, t)=Y^,i—i 5{r—fi(t)) at two different times At = t 2 — ii 
Q(At) = J d 3 ri d 3 r 2 p{r 1 ,0)p(r 2 , At) cos [qin-^)} 

N N 

= £E cos [^(o)-^(A*))]- (12) 

1=1 i=l 

Q(At) counts the number of "overlapping" particles in 
two configurations separated by a time interval At. Here, 
we do not apply a strict cutoff to define overlapping parti- 
cles, as was done in previous studiesj 2 ^* 2 ^^ but we follow 
Bcrthier™ and use the intermediate scattering function, 
cos [q(ri—r 2 )], where we again choose q — 1.7 A -1 . The 
fluctuations of this quantity are described by a general- 
ized susceptibility 

Xi (At) = ^[{Q*(At))-(Q(At)) 2 ], (13) 

which corresponds to the volume integral of a four-point 
density correlation function 

Xi(At) = ^ J d 3 r 1 d 3 r 2 d 3 r 3 d 3 r 4 cos [gift -7=3)] 

x cos [<f(r 3 -r 4 )] G 4 (fi, f 2 , f 3) r 4) At) (14) 
where /3 = fcsT and 

G4(ri,r2,r 3 ,r4,i) = (p(r 1 ,0)p(f 2 ,t)p(r 3 ,0)p(r 4 ,t)) 

-(p(r ll 0)p(r 2 ,t))(p(r 3 ,0)p(r i ,t)). 

The overlapping particles are comprised of localized par- 
ticles that have hardly moved and particles that have 
been replaced by a neighboring particle. Accordingly, Xi 
can be decomposed into self- (xss)i distinct- (xdd) and 
interference (xsd) parts, where the main contribution 
was found to result from the self parti^2& In this sec- 
tion, we are interested in spatial correlations of highly 
immobile particles. Therefore, we focus on the self part, 
which is given by 

Xss(At) = ^ [<Q|(Ai)) - (Qs(At)) 2 ] , (15) 
where 

N 

Qs(At) =J2cos[qXf f i(Q)-f i (At))}. (16) 

In Fig. we show xss(At) for the oxygen and sili- 
con atoms at several values of T. While the susceptibility 
is small at high T, it shows a pronounced maximum at 
low T. This strong increase of the peak height, xT> m ~ 
dicates a growing range of spatial correlations between 
immobile particles with decreasing 7 1 ^2^26^27 Unfortu- 
nately, a detailed analysis of the temperature dependence 
of xT i s hampered by a large scattering of the values 
(not shown). To characterize the times when the fluc- 
tuations in the number of localized particles are biggest, 
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FIG. 11: Pair correlation functions ffooM at T = 3030K. We 
compare results obtained for all oxygen atoms with those for 
the most mobile and the most immobile oxygen atoms in a 
time interval At=10ps~fs, respectively. 



we extract the peak times t±. From the results in Fig. [21 
it is evident that, for both atomic species, the spatial 
correlations between highly immobile particles are max- 
imum in the a-relaxation regime. In particular, (r) and 
ti show a very similar temperature dependence, where 
the ratio (r)/i4 amounts to a factor of about two. On 
the other hand, t± is about one order of magnitude longer 
than ts and, hence, SHD analyzed from the perspectives 
of highly immobile and highly mobile particles, respec- 
tively, is most pronounced at different times. All these 
findings are consistent with previous results for BLJ liq- 

24.25 .26.27.28 





average 


mobile 


zoo 


7.2 


7.5 


ZSiSi 


4.2 


4.4 


zosi(n = 2) 


0.99 


0.94 


zsio(n = 4) 


0.98 


0.94 



TABLE I: Coordination numbers characterizing the local 
structure in BKS silica at T = 3030K. 



nation number n is defined as the number of (3 atoms for 
which the distance from a given a atom is smaller than 
that corresponding to the first minimum of g a f}(r). It 
is evident that mobile oxygen and silicon atoms exhibit 
somewhat higher mean coordination numbers zoo and 
~zsiSi > respectively, suggesting that their positions are less 
favorable in terms of the potential energy. Further, the 
percentage of mobile oxygen atoms in the "ideal" coordi- 
nation with two silicon neighbors is smaller than for an 
average oxygen atom. Likewise, compared to an average 
particle, mobile silicon atoms are more often three- or 
fivefold coordinated to oxygen. However, the vast ma- 
jority of the mobile particles are still in their preferred 
local structure at the beginning of the particular time in- 
terval and, hence, defects of the local structure are not a 
necessary precondition of high particle mobility. On the 
other hand, the observed correlations imply that a some- 
what less ordered local environment facilitates particle 
dynamics. Thus, the relation between structure and dy- 
namics in BKS silica is non-trivial, as was also observed 
for BLJ mixtures 22 and the Dzugutov liquid. 29 



E. Relation between structure and dynamics 

Amorphous silica consists of a network of well-defined 
silicate tetrahedra. Hence, one may speculate that high 
particle mobility is found where the local structure is 
distorted. To study the relation between structure and 
dynamics, we identify the most mobile and the most im- 
mobile particles during a time interval At and character- 
ize their respective environments at the beginning of this 
period by means of the pair correlation functions (PCF), 
5a/3(^)j and the distributions of coordination numbers, 
z a 0{n) (a,(3£ {0,Si}). In Fig. ^2 we compare goo{r) 
for the mobile and immobile oxygen atoms, where At «tg 
and T = 3030 K. While the data for the immobile par- 
ticles are comparable to those for the ensemble average, 
the first and second neighbor peak are slightly broadened 
for the mobile oxygen atoms. This behavior is similar for 
all g a p{r). Thus, mobile particles during a time interval 
Ats exhibit a somewhat less ordered local environment at 
the beginning of this time window, at least as measured 
by g a i3{ r )i but the effects are weak. A variation of At in 
a meaningful range does not change these conclusions. 

Some of the results from our analysis of the distribu- 
tions z a p(n) are compiled in Table [IJ where the coordi- 



IV. DISCUSSION 

One goal of the present work is to compare dynamic 
heterogeneity in various model liquids. For this purpose, 
we show our results together with the corresponding lit- 
erature data in Table ITU To enable a quantitative com- 
parison, we consider results for T ~ Tmct in each case. 
Inspecting the data, it becomes clear that there is a corre- 
lation between the properties of dynamic heterogeneity 
at intermediate times and the features of the a relax- 
ation, e.g., the stretching parameter /3 characterizing the 
nonexponentiality of the decay of the incoherent inter- 
mediate scattering function. Specifically, the Dzugutov 
liquid shows the most pronounced dynamic heterogene- 
ity at intermediate times and the smallest value of (3, 
whereas systems like silica, with less heterogeneous dy- 
namics at intermediate times, exhibit a larger (3. In what 
follows, we discuss the findings for various model liquids 
in more detail. 

First, we see from Table ITT1 that the non-Gaussian pa- 
rameter 0L2 attains the largest values for the non-network 
forming liquids. Though ai was found to characterize as- 
pects of the heterogeneity of dynamics^ it does not pro- 
vide direct information about the spatial arrangement 
of mobile and immobile particles. The spatially het- 
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0.75 


0.55 


0.70 


0.50 


0.75 


0.83 


0.85 




1.6 
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1.3 


0.8 


Qm 

Jw 


15 (5) 




16(6.5) 


22 (5) 




7.9 (7) 


3.9(7) 


J W 






6.0 (6.5) 


11(5) 


4.2 (7) 


4.9 (7) 


2.8(7) 


r m 


2.2(5) 




1.9(6.5) 


2.4 (5) 




1.5(5) 


1.1(5) 


r 








0.70 (5) 




0.24(5) 


0.07(5) 


xT 


4.7 


18 













TABLE II: Quantities characterizing dynamic heterogeneity in various model liquids at TksTmct- We compare data for a 80:20 
binary Lennard- Jones mixture (BLJi, T= 1. 04 Ta/ct) a 50:50 binary Lennard- Jones mixture (BLJ2, T = Tmct),~^~ a 

polymer melt (PM, T= 1.02 Tmct) , 34 ' 35 ' 62 the Dzugutov liquid (DZ, T= 1.05 Tmct)^- and water (H2O, T= IMTmct)^^ 
with the present results for the oxygen (O) and silicon (Si) atoms in BKS silica at T = Tmct- The values in brackets denote 
the fractions of mobile particles, (f), used in the respective analysis (in percent). 



erogeneous nature of dynamics in the vicinity of Tmct 
can be studied by identifying clusters of mobile parti- 
clesi22i2L2ii2I On a qualitative level, the clusters in vari- 
ous model liquids, including BKS silica, show an analo- 
gous behavior. For example, the clusters strongly grow 
upon cooling, where, at each T, the mean cluster size is 
maximum at times tg in the late-/?/ early-a relaxation 
regime. On a quantitative level, several findings for sim- 
ple liquids22iSI and a polymer melt^*^ e.g., a power-law 
distribution of the cluster size, cannot be generalized to 
the case of silica. Moreover, there is a spectrum of max- 
imum mean cluster sizes, the higher and lower end of 
which are attained by the Dzugutov liquid and the net- 
work formers silica and water, respectively, see Tab. ITT1 
However, one has to consider that, with the same defini- 
tion, see section IlII Bl the sizes of the clusters depend on 
the properties of the local structure. Specifically, com- 
pared to simple liquids, the particles in network-forming 
liquids exhibit smaller coordination numbers and, hence, 
the probability that a mobile particle has a mobile neigh- 
bor is smaller. This means that, within the same defini- 
tion, the clusters are less likely to grow. At least to some 
extent, these effects are canceled out when the mean clus- 
ter size is normalized by the value expected from random 
statistics. However, at the present time, it is not clear 
whether Sn and Sw allow one to study SHD completely 
independent of the local structure. In any case, a com- 
parison of the results among the simple liquids, where the 
coordination numbers are similar, suggests a relation be- 
tween the mean cluster size and the stretching parameter 
P. 

A main difference between various model liquids is the 
relevance of string- like motion at T ~ Tmct- From the 
mean string length LJ} and the fraction f m of mobile par- 
ticles moving in non-trivial strings, it is evident that this 
type of motion is an important channel of relaxation in 
the Dzugutov liquid, but not in BKS silica, see Tab. ITU 
Moreover, the two atomic species in silica behave dif- 
ferently. On the one hand, some string-like motion is 
observed for the oxygen atoms. In agreement with pre- 
vious resultsj22i22ii^ the mean string length is maximum 
at intermediate times ti,~ts, where the peak height L 1 ^ 



increases with decreasing T. Hence, despite a limited rel- 
evance at TpzTmct, string- like motion may become im- 
portant for the structural relaxation of the oxygen atoms 
near T g , cf. Fig.0 On the other hand, the silicon atoms 
do not show string-like motion at intermediate times. In- 
stead, there are very few short strings at much later times 
ts "^tL~ (t) , implying that their nature is different from 
that observed for other model liquids ^^SiSk Therefore, 
we conclude that the "usual" string-like motion is ab- 
sent for the silicon atoms. When we also consider that 
the oxygen and silicon atoms typically exhibit two and 
four covalent bonds, respectively, our findings imply that 
string-like motion is suppressed by the presence of cova- 
lent bonds. 

Giovambattista et alM> found for water that the diffu- 
sion coefficient D is related to the mean cluster size via 
the AG relation where the cluster size is a measure of the 
mass of the CRR. In particular, they showed that a linear 
relation exists between log D and <5™— 1. For silica in the 
studied temperature range, we observe no such relation 
between D and either the mean cluster size or the mean 
string size and, hence, the findings for water cannot be 
generalized to the case of BKS silica. 

An analysis of the generalized susceptibility Xss(At) 
revealed spatial correlations of localized, i.e. highly im- 
mobile, particles in BKS silica. On a semi-quantitative 
level, our results are in good agreement with that for 
BLJ liquidsi 24 ' 25 ! 26 ! 27 More precisely, x ss (Ai) shows a 
maximum that strongly increases with decreasing T, in- 
dicating a growing length of spatial correlations between 
highly immobile particles. Moreover, the peak time ti is 
located in the cv-relaxation regime and closely follows the 
temperature dependence of the structural relaxation. A 
quantitative comparison of the present values of xT with 
literature data is hampered by the different units (re- 
duced/real) used in the simulations. However, the data 
for the 80:20 BLJ and 50:50 BLJ liquids again imply that 
the structural relaxation is more stretched for systems 
with pronounced SHD. 

Very recently, Garrahan and Chandler&S introduced a 
microscopic model of supercooled liquids, which is based 
on three central ideas: (i) Particle mobility is sparse and 
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dynamics are spatially heterogeneous at times interme- 
diate between ballistic and diffusive motion, (ii) Particle 
mobility is the result of dynamic facilitation, i.e., mobile 
particles assist their neighbors to become mobile, (iii) 
Mobility propagation carries a direction, the persistence 
length of which is larger for fragile than for strong glass 
formers. Our present and previous^ - results support the 
main ideas of the Garrahan-Chandler model on a quali- 
tative level. Specifically, we showed that SHD exists in 
BKS silica and, hence, this phenomenon is not limited to 
non-strong model liquids. Moreover, it was demonstrated 
in our previous work on BKS silica 40 that dynamic facili- 
tation is important, which is further corroborated by pre- 
liminary results for the Dzugutov liquid^ Finally, Gar- 
rahan and Chandler argue that the persistence length of 
the direction of mobility propagation manifests itself in 
the relevance of string-like motion so that the very lim- 
ited importance of this type of motion for BKS silica is 
consistent with their idea of a shorter persistence length 
of particle flow direction in strong liquids. 

We conclude that, at least for BKS silica, there is a 
subtle relation among dynamical processes on different 
time scales. Specifically, we suggest that the different 
properties of viscous liquids on the timescale of the struc- 
tural relaxation are not only a consequence of quantita- 
tive differences in SHD at intermediate times, but also 
result from differences in the spatiotemporal propagation 
of mobility. Consistently, we find that the evolution of 
SHD at intermediate times is insensitive to the crossover 
in the temperature dependence of the transport coeffi- 
cients, and vice versa. 

Finally, we comment on the relevance of our findings 
with respect to real silica. The BKS model is one of the 
simplest models of silica, and neglects atomistic details 
such as three-body interactions and charge transfer pro- 
cesses, which are included in newer silica potentials ' 4 ^ 
Both of these features have been demonstrated to be im- 
portant for, e.g., structural transformations between sil- 
ica crystal polymorphs^, Despite its simplicity, the BKS 
model reproduces many bulk dynamic and thermody- 
namic features of real silica. However, local dynamics 
involving correlated or string- like motion in the melt may 
well be sensitive to the re-distribution of charge during 
bond breakage and to higher order terms in the inter- 
action potential. Investigation of SHD in more realistic 
silica models would determine this, and are underway. 

V. SUMMARY 

We showed that dynamics in BKS silica are spatially 
heterogeneous and, hence, the structural relaxation in 
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this model of a strong liquid cannot be described as a 
statistical bond-breaking process, as may have been ex- 
pected from the Arrhenius-like temperature dependence 
of the transport coefficients. Specifically, we demon- 
strated that there are spatial correlations between mobile 
and immobile particles, respectively, the extent of which 
grows strongly upon cooling. In particular, the growth 
of dynamically correlated regions continues in the Arrhe- 
nius temperature regime. While the spatial correlations 
between highly mobile particles are maximum in the late- 
(5/ early-a relaxation regime of the MCT, those between 
highly immobile particles are biggest at later times close 
to the time constant of the structural relaxation. On a 
qualitative level, all these findings for SHD in BKS sil- 
ica resemble that for non-strong model liquids. Consis- 
tently, we found that defects of the local network struc- 
ture facilitate dynamics to some extent, but they are not 
a necessary precondition for high particle mobility. On a 
quantitative level, a detailed comparison with literature 
data showed that measures of dynamic heterogeneity at 
intermediate times, e.g., the mean size of clusters of mo- 
bile particles, exhibit a broad spectrum of values at the 
lower end of which the present results for BKS silica are 
found. Moreover, such comparison revealed that the a- 
relaxation is more nonexponential for liquids with pro- 
nounced dynamic heterogeneity at intermediates times. 
However, we found no evidence for a straightforward re- 
lation between the properties of dynamics on different 
time scales so that the correlations are subtle. In partic- 
ular, the "fragile-to-strong crossover" for the transport 
coefficients is not accompanied by a substantial change 
in the temperature dependence of SHD at intermediate 
times. On the other hand, our results are qualitatively 
consistent with a microscopic model of viscous liquids put 
forward by Garrahan and Chandler. 8,9 Following their 
ideas, we suggest that the spatiotemporal characteristics 
of mobility propagation play an important role in the 
differences between fragile and strong liquids. 



Acknowledgments 



We thank Y. Gebremichael and M. Bergroth for stimu- 
lating discussions. M. V. gratefully acknowledges funding 
by the Deutsche Forschungsgemeinschaft (DFG) through 
the Emmy Noether-Programm. 



many; Electronic address: mivogel@uni-muenster.de 
Electronic address: sglotzer@umich.edu 



13 



1 K. Schmidt-Rohr and H. W. Spiess, Phys. Rev. Lett. 66, 
3020 (1991) 

2 M. T. Cicerone, F. R. Blackburn and M. D. Ediger, Macro- 
molecules 28, 8224 (1995) 

3 R. Bohmer, R. V. Chamberlin, G. Diezemann, B. Geil, A. 
Heuer, G. Hinze, S. C. Kuebler, R. Richert, B. Schiener, 
H. Sillescu, H. W. Spiess, U. Tracht and M. Wilhelm, J. 
Non-Cryst. Solids 235-237, 1 (1998) 

4 H. Sillescu, J. Non-Cryst. Solids 243, 81 (1999) 

5 M. Ediger, Ann. Rev. Phys. Chem. 51, 99 (2000) 

6 L. A. Deschenes and D. A. Vanden Bout, Science 292, 255 
(2001) 

7 G. Adam and J. H. Gibbs, J. Chem. Phys. 43, 139 (1965) 

8 J. P. Garrahan and D. Chandler, Phys. Rev. Lett. 89, 
035704 (2002) 

9 J. P. Garrahan and D. Chandler, Proc. Nat. Acad. Soc. 
100, 9710 (2003) 

10 U. Tracht, M. Wilhelm, A. Heuer, H. Feng, K. Schmidt- 
Rohr and H. W. Spiess, Phys. Rev. Lett. 81, 2727 (1998) 

11 S. A. Reinsberg, X. H. Qiu, M. Wilhelm, H. W. Spiess and 
M. D. Ediger, J. Chem. Phys. 114, 7299 (2001) 

12 X. H. Qiu and M. D. Ediger, J. Phys. Chem. B 107, 459 
(2003) 

13 W. Gotze, L. Sjogren, Rep. Prog. Phys. 55, 241 (1992) 

14 D. N. Perera and P. Harrowell, Phys. Rev. Lett. 81, 120 
(1998) 

15 Y. Hiwatari and T. Muranaka, J. Non-Cryst. Solids 235- 
237, 19 (1998) 

16 A. Onuki and R. Yamamoto, J. Non-Cryst. Solids 235-237, 
34 (1998) 

17 B. Doliwa and A. Heuer, Phys. Rev. Lett. 80, 4915 (1998) 

18 G. Wahnstrom, Phys. Rev. A 44, 3752 (1991) 

19 G. Johnson, A. I. Mel'cuk, H. Gould, W. Klein and R. D. 
Mountain, Phys. Rev. E 57, 5707 (1998) 

20 W. Kob, C. Donati, S. J. Plimpton, P. H. Poole and S. C. 
Glotzer, Phys. Rev. Lett. 79, 2827 (1997) 

21 C. Donati, J. F. Douglas, W. Kob, S. J. Plimpton, P. H. 
Poole and S. C. Glotzer, Phys. Rev. Lett. 80, 2338 (1998) 

22 C. Donati, S. C. Glotzer, P. H. Poole, W. Kob, S. J. Plimp- 
ton, Phys. Rev. E 60, 3107 (1999) 

23 C. Donati, S. C. Glotzer and P. H. Poole, Phys. Rev. Lett. 
82, 5064 (1999) 

24 S. C. Glotzer, V. N. Novikov and T. B. Schr0der, J. Chem. 
Phys. 112, 509 (2000) 

25 N. Lacevic, F. W. Starr, T. B. Schr0der, V. N. Novikov 
and S. C. Glotzer, Phys. Rev. E 66, 030101 (2002) 

26 N. Lacevic, F. W. Starr, T. B. Schr0der and S. C. Glotzer, 
J. Chem. Phys. 119, 7372 (2003) 

27 L. Berthier, Phys. Rev. E 69, 020201 (2004) 

28 C. Donati, S. Franz, S. C. Glotzer and G. Parisi, J. Non- 
Cryst. Solids 307, 215 (2002) 

29 M. Dzugutov, S. I. Simdyankin and F. H. M. Zetterling, 
Phys. Rev. Lett. 89, 195701 (2002) 

30 Y. Gebremichael, M. Vogel and S. C. Glotzer, J. Chem. 
Phys. 119, 4415 (2004) 

31 Y. Gebremichael, M. Vogel and S. C. Glotzer, Molecular 
Simulation 30, 281 (2004) 

32 A. Heuer and K. Okun, J. Chem. Phys. 106, 6176 (1997) 

33 C. Bennemann, C. Donati, J. Baschnagel and S. C. 
Glotzer, Nature 399, 246 (1999) 



34 Y. Gebremichael, T. B. Schr0der, F. W. Starr and S. C. 
Glotzer, Phys. Rev. E 64, 051503 (2001) 

35 M. Aichele, Y. Gebremichael, F. W. Starr, J. Baschnagel 
and S. C. Glotzer, J. Chem. Phys. 119, 5290 (2003) 

36 J. Qian, R. Hentschke and A. Heuer, J. Chem. Phys. 110, 
4514 (1999) 

37 N. Giovambattista, S. V. Buldyrev, F. W. Starr and H. E. 
Stanley, Phys. Rev. Lett. 90, 085506 (2003) 

38 S. C. Glotzer, J. Non-Cryst. Solids 274, 342 (2000) 

39 C. Dasgupta, A. V. Indrani, S. Ramaswamy and M. K. 
Phani, Europhys. Lett. 15, 307 (1991) 

40 M. Vogel and S. C . Glotzer, Phys. Rev. Lett, (in press); 
|cond-mat/0402427| 

41 C. A. Angell, J. Non-Cryst. Solids 131-133, 13 (1991) 

42 B. W. H. van Beest, G. J. Kramer, R. A. van Santen, Phys. 
Rev. Lett. 64, 1955 (1990) 

43 J.-L. Barrat, J. Badro and P. Gillet, Molecular Simulation 
20, 17 (1997) 

44 J. Horbach and W. Kob, Phys. Rev. B 60, 3169 (1999) 

45 J. Horbach and W. Kob, Phys. Rev. E 64, 041503 (2001) 

46 W. Kob, J. Phys.: Condens. Matter 11, R85 (1999) 

47 K. Binder, J. Non-Cryst. Solids 274, 332 (2000) 

48 I. Saika-Volvod, P. H. Poole and F. Sciortino, Nature 412, 
514 (2001) 

49 K. Vollmayr, W. Kob and K. Binder, Phys. Rev. B 54, 
15808 (1996) 

50 J. Horbach, W. Kob, K. Binder and C. A. Angell, Phys. 
Rev. E 54, 5897 (1996) 

51 J. Horbach, W. Kob and K. Binder, Eur. Phys. J. B 19, 
531 (2001) 

52 K.-U. Hess, D. B. Dingwell and E. Rossler, Chem. Geol. 
128, 155 (1996); E. Rossler, K.-U. Hess and V. N. Novikov, 
J. Non-Cryst. Solids 223, 207 (1998) 

53 M. Bergroth, M. Vogel and S. C. Glotzer (unpublished 
results) 

54 M. P. Allen and D. J. Tildesley, Computer Simulation of 
Liquids (Oxford University Press, New York, 1990) 

55 J. C. Mikkelsen, Appl. Phys. Lett. 45, 1187 (1984) 

56 G. Brebec, R. Seguin, C. Sella, J. Bevenot, J. C. Martin, 
Acta Metall. 28, 327 (1980) 

57 F. W. Starr, F. Sciortino and H. E. Stanley, Phys. Rev. E 
60, 6757 (1999) 

58 S. Sastry, Nature 409, 164 (2001) 

59 A. Scala, F. W. Starr, E. La Nave, F. Sciortino and H. E. 
Stanley, Nature 406, 166 (2000) 

60 Y. Gebremichael, PhD thesis, University of Maryland, 
USA (2004); Y. Gebremichael, M. Vogel, F. W. Starr and 
S. C. Glotzer (in preparation) 

61 T. B. Schr0der, S. Sastry, J. C. Dyre and S. C. Glotzer, 
J. Chem. Phys. 112, 9834 (2000); T. B. Schroder, 
cond-mat/0005127 

62 C. Bennemann, W. Paul, K. Binder and B. Diinweg, Phys. 
Rev. E 57, 843 (1998) 

63 F. W. Starr, S. Harrington, F. Sciortino and H. E. Stanley, 
Phys. Rev. Lett. 82, 3629 (1999) 

64 L. Duffrene and J. Kieffer, J. Phys. and Chem. of Solids 
59, 1025 (1998) 

65 L. Huang and J. Kieffer, J. Chem. Phys. 118, 1487 (2003) 



